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Abstract 

We show that the jet structure of DM annihilation or decay prod- 
ucts enhances the d production rate by orders of magnitude com- 
pared to the previous computations done assuming a spherically 
symmetric coalescence model. In particular, in the limit of heavy 
DM, M ^ TTT-p, we get a constant rather than suppressed 
d production rate. Therefore, a detectable d signal is compatible 
with the lack of an excess in the p PAMELA flux. Most impor- 
tantly, cosmic d searches become sensitive to the annihilations 
or decays of heavy DM, suggesting to extend the experimental d 
searches above the 0{1) GeV scale. 
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1 Introduction 



The cosmological DM abundance is naturally produced in thermal freeze-out if Dark 
Matter (DM) has weak interactions and a TeV-scale mass M, that in appropriate models 
can be lowered down to the weak scale, 100 GeV. This scenario can be tested searching for 
DM annihilation (or decay) products in cosmic rays. In view of astrophysical backgrounds, 
a particularly sensitive signal is an excess in cosmic-ray anti-particles: positrons, anti- 
protons p and anti-deuterium d. According to the coalescence prescription [1], ad is 
formed when DM produces a p and a n with momentum difference below Pq ~ 160 MeV. 
The standard formula for the d spectrum, obtained under the assumption of spherical 
symmetry of the events, in terms of the anti-nucleon {p plus n) energy spectrum per 
annihilation, dNj^/dT, is [2, 3, 4, 6, 7] 

dN^_ pI fldNA' 
dTd 3kjmp \2 dT ) T=Ta/2 

where the kinetic energies T = E — m are Tp = Tn = Trf/2, rrip = rrin = rnd/2 and 
/cj = \Jt^ + 2mdTj. Eq. (1) implies a d yield suppressed by 1/M^ for large M. This 
result is qualitatively wrong. Increasing M just increases the boost of the primary DM 
annihilation products, giving rise, due to Lorentz symmetry, to an essentially constant d 
production rate with energy roughly proportional to M. The reason for this fundamental 
discrepancy is caused by the fact that the spherical approximation misses the jet structure 
of the DM annihilation products. 

In this letter we show that for M ^ nip the angular proximity of the produced 
p, n enhances the d yield, possibly by orders of magnitude. We critically compare the 
standard spherical approximation results with our Monte Carlo approach to d produc- 
tion, presenting the d energy spectra for the various DM annihilation or decay channels 
into W~^W~ , ZZ, qq, bb, ti, hh and comment on the astrophysical d background produced 
mostly in cosmic ray pp and pp collisions. We propagate d in the Milky Way, studying the 
phenomenology and the prospects for DM produced d searches at AMS-2, in the light of 
the PAMELA p observations. We find that the d signal below 1 GeV is strongly enhanced 
increasing the chances of d detection at AMS-2 even for the standard thermal DM anni- 
hilation cross section. This result is consistent with the lack of p/p excess in PAMELA. 
Due to the qualitatively different large M behavior of the production rate, our result 
drastically enhance the d production at high energies. Therefore the cosmic ray d flux 
produced in heavy DM annihilations or decays exceeds the estimated background, and 
AMS-2 and future d experiments become sensitive to DM if they extend their sensitivity 
to d above 1 GeV. 
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2 Spherical-cow vs Monte Carlo 



d is formed via pn —>■ d'j, that has a large cross-section due to the small binding energy 
of d, which therefore has a spatially extended wave-function or equivalently a strongly 
peaked wavefunction iplAk) = {d\pn) in momentum space. Here d has momentum kd and 
energy Ed and Ak is the relativistically invariant relative momentum between p (with 
momentum kp and energy Ep) and n (with momentum fc„ and energy En)'- 

AP = \kp - - {Ep - EnY + (m„ - nipY, (2) 

where we can neglect rrip — rrin = 1.29 MeV. The amplitude for d production in DM 
annihilations, DM DM —* d, can be computed as 

((J|DM DM) = J2{d\pn){pn\DM DM), (3) 



p,n 

|2 



giving the 'coalescence approximation' [1]: the probability |((i|pn)| that a p and a n 
coalesce to form a c? is approximated as a narrow step function 6(AA; — po), that drops 
from unity to zero if Ak is larger than po- Here pq is a constant (to be extracted from 
data later) that can be estimated as po ~ y/rndBd ~ 60 MeV assuming that d production 
happens until the relative p, n kinetic energy is smaller than the deuteron binding energy 
Bd = 2.2 MeV. The total d yield therefore is 

Nd = J dNp dNM^e-pl) = I d%d'kn^^^^@{Ae-pl). (4) 

In the non-relativistic limit kp^n ^p,n and for small po the region that satisfies AA; < pq 
at fixed A;„ is a sphere in kp centered on k^ with radius Pq and volume Airp^/?). In general 
the sphere gets dilatated along the direction kp ^ /c„ by a relativistic Lorentz factor 
7p ~ 7n ~ Id- Multiplying eq. (4) times 1 = / d^kd5{kd — kp — k^) we finally get, in the 
limit of small Pq -C M/'^p^ni the d momentum distribution: 

dN^_lA7Tp^ dNpdNn 
^U^kd 8 3 ^""^"d^kpd^kn' 

where kp = kn = kd/'2. Eq. (5) is relativistically invariant as it contains the usual 
relativistic phase space d^k/2E = d'^k 5{E'^ — k"^ — th?). 

The spherical approximation 

Previous computations proceed assuming spherical symmetry, d^k = Arrk'^ dk, and uncor- 

related p, n distributions: 

dNpdNn dNp dNn , E dN 1 dN 

implymg = t-t—^t^- (6) 



d^kpd^kn d^kp d^kn md^k Airkm dE 



^The extra factor of 8 with respect to the equation used in papers [2]- [7] comes from d^kd — Sd^kp^n- 
In the final result this difference gets compensated by a value of po twice larger than the one adopted in 
those papers. In our Monte Carlo computation of the coalescence condition Afc < pq it is important that 
we fix factors of 2 so that our po really is the radius of the coalescence sphere. 
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Writing the result in terms of the adimensional Xi = Ti/M (so that < Xd,p,n < 1 and 
Xp = Xn = Xd/^) one gets eq. (1) i.e. 



which is exphcitly suppressed by for large DM mass M. 

This is qualitatively wrong. Consider for example the DM DM W'^W~ annihilation 
mode. Increasing M increases the boost of each W, and thereby the boost of the anti- 
deuterons from W decay, but the d number stays fixed. Neglecting QED final state 
radiation (FSR), for M ^ M\y one should get a constant, M-independent function for 
dNd/dxd- Obviously the problem is in the 'spherical cow' approximation [8]. Due to the 
boost the events are highly non spherical and SM particles are concentrated in two 
back-to-back jets, enhancing the probability of having pn pairs with small momentum 
difference Ak < pq. A similar argument applies to DM annihilations or decays into 
colored particles, such as qq. Hadronization leads to QCD jets, rather than to spherical 
events. Thereby the spherical approximation can grossly underestimate the d production. 

Going to less relevant aspects that control order one factors, the analytic spherical 
approximation can also over-estimate the d yield, by neglecting anti-correlations between 
n and p or the fact that no d is obtained if only one anti-nucleon is present per event. As 
an example, we consider again the W^W~ mode: within the spherical approximation a d 
can form coalescing a p from W~ with a n from , but this process is highly suppressed 
because the and W~ go back to back. 

The Monte Carlo approach 

In order to take into account the jet structure of the events and the correlations between 
the p,n momenta we compute the d spectrum by searching event-by-event for the n,p 
pair(s) which have relativistically invariant momentum difference Ak smaller than pq. 
We verified that the spherical uncorrelated approximation of eq. (7) is reproduced if we 
first merge many events, and later coalesce p with n without imposing that they come 
from the same event. 

Various experiments extracted compatible values of po from data about d production 
in hadronic and e~^e~ colhsions. Presumably these studies adopted the 'spherical cow' 
approximation rather than performing a Monte Carlo computation. Giving the relatively 
low energies involved this should not make a large difference; anyhow we here prefer to 
directly extract po from the ALEPH data [9]: one hadronic Z decay at rest gives rise 
to (5.9 ± 1.9) 10~^ d in the momentum range 0.62 GeV < kd < 1.03 GeV and angular 
range |cos^| < 0.95. According to our Monte Carlo computation, this translates into 
Po = 162 ± 17 MeV. Should po have a value different from the po = 160 MeV adopted here 
for both the DM signal and the astrophysical background (as computed in [2, 4]) the d 
energy spectra get rescaled roughly by an overall (po/160MeV)^ factor. 




(7) 
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Figure 1: Total d yield per DM annihilation as a function of the DM mass. The thick 
upper lines present the Monte Carlo result, and the lower thin lines are the spherical 
approximation. The annihilation modes are into W^W~ (red), ZZ (blue), hh (green) tt 
(black dot-dashed), bb (blue dashed), light quarks qq (red dashed). 

We performed a Monte Carlo study by generating a huge number of events (up to 10'^ 
per DM DM annihilation, and up to 10^ events when studying pp and pp collisions) with 
Pythia 8 [10], directly implementing the condition Afc < po = 160 MeV for d production. 
Such computing-power demanding results have been obtained using the EU Baltic Grid 
facilities [11]. 

Fig. 1 shows the total number of d produced per DM annihilation as function of the DM 
mass for various annihilation modes, comparing our Monte Carlo result with the spherical 
approximation, which can under-estimate the d yield by various orders of magnitude. The 
same po = 160 MeV is assumed in both cases. 

Fig. 2 shows our Monte Carlo results for the d spectra computed for three values 
of the annihilating DM mass M. The same spectra also hold for decaying DM, after 
replacing Mann Mdec/2. As we expected, the result has only a minor dependence on M 
and is thereby qualitatively different from the 'spherical-cow' approximation that would 
give a 1/M^ suppression. There are three classes of qualitatively different cases: DM 
annihilations i) into W,Z,h (we assume a Higgs mass nih = 120 GeV); ii) into quarks 
q, b, t or iii) into leptons. The latter case gives no d. To compare the former two cases 
that give d, we focus on i) DM DM — > W^W~ and ii) DM DM qq, and show the d 
spectra in fig. 3a and b, respectively, for various values of the DM mass M. In the W^W~ 
case the d spectrum only mildly depends on the DM mass. Neglecting FSR, all d should 
have X > rrid/Mw = 0.05; the small d flux at smaller x is due to electroweak FSR. In the 
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Figure 2: Monte Carlo results for the d spectra x ■ dN/dx produced per DM annihilation 
into ZZ, hh tt, bb, qq. We assumed the DM mass M = 200 GeV (1 TeV ; /5 TeV/ 

in the left (middle) [right] panel. The notation is the same as in fig. 1. 



qq case the d spectrum at smaller x increases with M rather than being suppressed as 
This is due to QCD FSR that roughly scales as ln(M/mp). 
The Monte Carlo results differ both qualitatively and quantitatively from the previous 
studies of the d spectra in DM annihilations or decays. To draw conclusions about the 
detectability of the signal, we also need to study possible changes in the astrophysical d 
background, mainly generated by collisions of cosmic-ray p with energy Ep on p at rest. 
In view of the kinematical threshold for d production {Ep > 30 GeV) and of the energy 
spectrum of cosmic protons (roughly proportional to E'^), d production is dominated 
by Ep ~ 60 GeV, in the range of validity of the parton model in Pythia. Our semi- 
quantitative results for the d background suggest a reasonable agreement with the spectra 
of [2,4]. This is an expected result because the center of mass energy in cosmic pp collisions 
is small and, in this case, the uncorrelated spherical approximation is expected to work 
reasonably well. However this issue needs to be precisely investigated. 

Some remarks are in order. First, we computed p,n allowing all other hadrons to 
decay despite that the life-time of some strange baryons, such as the S = uss, is longer 
than the size of deuterium. This effect should already have been taken into account when 
extracting the value of Pq from high-energy experimental data from its definition of eq. (5). 
Second, DM in general annihilates into various primary channels k. According to eq. (1) 
one should sum their contributions to the p,n spectra rather than to the d spectrum, 
getting {J2kdN^'^ydx)'^ ^ Y.kidN'^^^ dxY . Our Monte Carlo result instead amounts to 
sum incoherently over all primary annihilation channels k as well as all secondary and 
tertiary contributions in the decay chain. 
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d kinetic energy x = T/M d kinetic energy x = T/M 

Figure 3: Monte Carlo results for the d spectra per DM annihilation, x ■ dN/dx for 
Po = 160 MeV. We consider DM masses M = 0.1, 1, lOTeV (black, blue, red continuous 
curves) and M = 0.3,3,30TeV (black, blue, red dashed curves) and the indicated DM 
annihilation modes. 



3 Cosmological fluxes 

To compute the dflux in the solar system we consider four possible Milky Way DM density 
profiles p(r) [12]: 

(1 + rQ/r^)/(l + r^/r^) isothermal, = 5kpc 

p(j;l ^\ {rQ/r){l + rQ/rs)y{l + r/rs)^ NFW, = 20kpc 

Pq I (rQ/r)i-i6(l + r0/r,)V(l + r/r,)i-84 Moore, r, = 30 kpc 

exp(— 2[(r/r5)° — {rQ/rs)°']/a) Einasto, = 20kpc, a = 0.17, 

keeping fixed the local DM density p(r = Vq) = p© = 0.3 GeV/ cm'^. Concerning diffusion 
of charged d in the galaxy, we approximate the diffusion region as a cylinder with height 2L 
centered on the galactic plane, a constant diffusion coefficient K = KqE^ and a constant 
convective wind directed outward perpendicularly to the galactic plane. We consider the 
min, med, max propagation models [13] for p, d, which are characterized by the following 
astrophysical parameters, 

Kq in kpc^/Myr L in kpc Konv in km/s 
e; n nm i i ^ c; 

(9) 



Finally, one must take into account annihilations of d on interstellar protons and Helium 
in the galactic plane (with a thickness of h = 0.1 kpc -C L) with rate Fann [2]. The 
solution to the diffusion equation for the energy spectrum of the d number density, /, 

- K{T) ■ V V + £ (sign(^) / Konv) = Q-2h 6{z) F,,,/, (10) 
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Figure 4: The p (left) and d (right) astrophysical functions , R(T) and Rd(T) of eq. (11), 
computed under different assumptions. In both cases, the dashed (solid) [dotted] bands 
assumes the min (med) [max] propagation configurations. Each band contains 4 lines, 
that correspond to the isothermal (red lower lines), NFW (blue middle lines), Einasto 
(magenta) and Moore (green upper lines) DM density profiles. 

acquires a simple factorized form in the "no-tertiaries" approximation that we adopt. The 
d flux in the galactic medium around the solar system can be written as 



fully analogous to the solution for the p flux in [14]. The function dNj/dT contains the 
particle physics input and was computed in the previous section. The function Rd{T) 
encodes the Milky Way astrophysics and is plotted in fig. 4b for various halo and prop- 
agation models. It roughly is some average containment time in the diffusion cylinder, 
and we verified that d generated outside it provide a negligible extra contribution even 
in the min scenario, where most DM annihilations occur outside the diffusion cylinder: 
the probability of re-entering is sizable, but the probability of diffusing up to the solar 
system is small. Going from DM annihilations to DM decays with life-time r one just 
needs to replace in eq. (11) {(Tv)pq/2M'^ with Pq/Mt; we do not plot the corresponding 
Rd{T) functions for DM decay as they essentially coincide with the R^ function for DM 
annihilations and the isothermal profile plotted in fig. 4b. Indeed, for all the considered 
DM profiles, DM decays close to the galactic center do not significantly contribute to 
the d flux at Earth, as for DM annihilations with the quasi-constant isothermal density 
profile. 

We notice that although RdiT) is significantly uncertain (especially below a few GeV), 




(11) 
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Figure 5: Assuming the av = 3 ■ 10~^^ cm^/sec suggested by cosmology, the NFW profile, 
MED propagation and DM masses M = {0.1, 1, 10} TeV, we compare the d flux obtained 
from the full computation (continuous lines) with the one from the spherical approxima- 
tion. The dotted line is the expected astrophysical background. 



the ratio with the corresponding astrophysical function R{T) for p is essentially fixed, so 
that the non-observation of a DM p excess puts robust bounds on the possible DM d flux. 
Indeed the p flux has been observed below 100 GeV by PAMELA [15] and agrees with 
astrophysical expectations, which are believed to have an uncertainty of about ±20% [4, 5] . 

Finally, we take into account the solar modulation effect, relevant only for non- 
relativistic d: the solar wind decreases the kinetic energy T of charged cosmic rays such 
that the energy spectrum d^^^i^/dT^ of J that reach the Earth with energy T® is approx- 
imatively related to their energy spectrum in the interstellar medium, d^^/dT, as [18] 

c^$,-e _ 2m,Te + T| d^j 

iT^- 2m,T + T^ dT' T-T^ + ec^,. (12) 

The so called Fisk potential (pp parameterizes in this effective formalism the kinetic energy 
loss. We assume (pp = 0.5 GV i.e. ecpp = 0.5 GeV. 



4 Results 

Fig. 5 compares our Monte Carlo results for the d flux with the spherical approxima- 
tion.^ The shading indicates the enhancement. We here assumed the NFW profile, MED 
propagation, and the DM annihilation cross section av = cifcosmo = 3 10~^^ cm'^/sec that 
reproduces the cosmological DM abundance via thermal freeze-out. 

^ Numerical results in some previous computations apparently included spurious factors of 2 related 
to dNp/dx{a,ftei neutron decay) w 2d-A'p /da; (before neutron decay) (this explains a discrepancy with [4]) 
and to GeV/nuc = GeV/2 (that affects the measure dT in d^^/dT; when comparing our plots with ones 
in previous papers, notice that we plot d^^/dlnT rather than d^g/dT). 
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Figure 6: The p/p ratio, for the same DM models described in the caption of fig. 7, 
showing that they are compatible with the PAMELA p data (red points). Shading indicates 
the expected astrophysical background. 




d kinetic energy T in GeV/n d kinetic energy T in GeV/n 



Figure 7: Upper row: our result for the d flux at Earth. Lower row: previous re- 
sults for the d flux computed in the spherical approximation. We consider DM masses 
M = 0.1,1, 10 TeV (black, blue, red continuous curves) and M = 0.3,3,30TeV (black, 
blue, red dashed curves), DM annihilations into W~^W~ (left) and qq (right) with 
av = 3- 10-2^ cmVsec x max(l, M/300 GeV)^ the NEW DM profile, MED propagation, 
solar modulation (pp = 0.5 GV, po = 160 MeV. Shading indicates the expected astrophysi- 
cal d background. 



Rather than relying on theoretical assumptions, in order to explore the maximal d flux 
from DM compatible with present data, we assume 

av = max(l, M/300 GeV)^ ■ avcosmo • (13) 

Indeed, fig. 6 shows that these assumptions give a p fiux compatible with (and comparable 
to) the PAMELA p/p data. As discussed in the previous section, the p/d ratio is negligibly 
affected by astrophysical uncertainties. Furthermore, the assumed cross section is about 
one order of magnitude below what is needed to explain the PAMELA [19] e"*" excess and 
is compatible with the bounds from galactic 7 and u observations [20] as well as with the 
diffused 7-ray constraints [21]. 

Then, the upper row of fig. 7 shows our Monte Carlo results for the d fiux, T ■ d^^i/ dT, 
while the lower row shows the corresponding much lower d fiux obtained in the spherical 
approximation. The caption describes all various assumptions. 

Comparison of these two results shows that the signal is enhanced in our Monte Carlo 
d computation: almost an order of magnitude for small d kinetic energies (T ~ 1 GeV) or 
lighter DM (M ~ 100 GeV) and orders of magnitude at higher energies or for heavier DM. 
Such enhancement does not depend on the assumed value of av Before our calculation it 
was believed that only the sub-GeV energy region is suitable for searches of a DM-induced 
d signal. Our result implies that heavy DM, as suggested by PAMELA and FERMI data, 
also induces detectable d signal at high energies. Therefore our result has important 
implications on the strategy of DM searches using the d signal. 

The red line in fig. 7a roughly shows the Minimal Dark Matter [14] prediction: DM 
with M ^ 10 TeV that makes Sommerfeld-enhanced annihilations into W~^W~ , giving rise 
to p and consequently to d at energies (per nucleon) above rUpM/Mw, not yet explored 
by PAMELA. 

The PAMELA [19], FERMI [22] and HESS [23] excesses suggest a DM inter- 
pretation in terms of multi-TeV DM that annihilates dominantly into leptons with a 
Sommerfeld-enhanced cross section [17, 24]. An interesting class of models with these 
properties is obtained by assuming that DM annihilates into a new vector with mass 
m < 2mp, that subsequently can only decay into the lighter e,/i, vr [25]. We notice that 
this condition is not strictly necessary neither for the Sommerfeld enhancement nor for 
compatibility with PAMELA p data: indeed if m > 2mp one would obtain p with energy 
larger than mpM/m, where M/m is the boost factor of the new vector. This boost is 
large enough not to give an unseen p excess below 100 GeV (the energy range explored 
by PAMELA so far) even if m is several tens of GeV, as in [17]. Similarly to the Minimal 
Dark Matter case, these models would give a fiux of d above 100 GeV. Our enhanced d 
signal should also be used to re-evaluate prospects of discovering supersymmetric Dark 
Matter candidates, which often annihilate into the W^W~ or hh modes we considered. 



11 



5 Conclusions 



We computed the d flux at Earth produced by DM annihilations or decays in the Milky 
Way using an event-by-event Monte Carlo technique run on the Grid, improving on 
previous computations that assumed spherically symmetric events and obtained a 
suppression of the d yield for heavy DM masses M. Due to the jet structure of high 
energy events implied by relativity no such suppression is present, and the d signal is 
strongly enhanced: by orders of magnitude for d energies above 10 GeV or DM masses 
above 1 TeV, as illustrated in fig. 7. The d astrophysical background seems not to be 
significantly affected, being dominantly generated by low-energy cosmic ray collisions. 
While the p and d fluxes suffer from significant astrophysical uncertainties, their ratio is 
robustly predicted. Thereby the non-observation of a p excess in PAMELA data implies 
an upper bound on the d DM flux. In the light of our enhanced d fluxes, we find that a d 
DM signal is still possible. For example, heavy DM models [14, 25] that can account for 
the PAMELA excess can lead to p and d excesses above 100 GeV/nucleon. 

Most importantly, our result implies that the experiments searching for cosmic ray d 
become sensitive to M > TeV mass DM, provided that the DM annihilation cross section 
is larger than what naively suggested by thermal freeze-out. Therefore it is important to 
extend future searches for d above the GeV energy range. For the moment, the AMS-2 
experiment is expected to achieve a very energy-dependent efficiency to d detection, so 
that AMS-2 would have a sensitivity to a flux down to 5 10~'''/(m^ sec srGeV/nuc) 
in the energy ranges 0.2GeV/nuc < T < IGeV/nuc (where time-of-flight is enough to 
discriminate d from p) and 2GeV/nuc < T < 4GeV/nuc (where the magnetic spectrom- 
eter is needed) [26]. According to previous d DM computations based on the spherically 
symmetric approximation, only the lower energy range was promising for DM searches. 
We have shown that the DM signal can manifest itself also at higher energies, where it is 
less affected by astrophysical uncertainties. 
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